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ABSTRACT 

Samples returned from comet 81P/Wild 2 by the 5'tarrfust mission confirm that substantial quantities 
of crystalline silicates were incorporated into the comet at the time of its formation. We investigate the 
constraints that this observation places upon protoplanetary disk physics, under the assumption that 
outward transport of particles processed at high temperatures occurs via a combination of advection 
and turbulent diffusion in an evolving disk. We also look for possible constraints on the formation 
locations of such particles. Our results are based upon one-dimensional disk models that evolve with 
time under the action of viscosity and photoevaporative mass loss, and track solid transport using an 
ensemble of individual particle trajectories. We find that two broad classes of disk model are consistent 
with the Stardust findings. One class of models features a high particle diffusivity (a Schmidt number 
Sc < 1), which suffices to diffuse particles up to 20 /im in size outward against the mean gas flow. For 
Sc > 1 such models arc unlikely to be viable, and significant outward transport appears to require that 
the particles of interest settle into a midplane layer that experiences an outward gas fiow. In cither class 
of models, the mass of inner disk material that reaches the outer disk is a strong function of the initial 
compactness of the disk. Hence, models of grain transport within steady-state disks underestimate 
the efficiency of outward transport. Neither model results in sustained outward transport of very 
large particles exceeding a millimeter in size. We show that in most circumstances, the transport 
efficiency falls off rapidly with time. Hence, high-temperature material must be rapidly incorporated 
into icy bodies to avoid fallback to small radii. We suggest that significant radial transport may only 
occur during the initial phase of rapid disk evolution. It may also vary substantially between disks 
depending upon their initial mass distributions. We discuss how our model may inform recent Spitzer 
observations of crystalline silicates in T Tauri star-disk systems. 

Subject headings: accretion, accretion disks — comets: individual (81P/Wild 2) — planets and satel- 
lites: formation — protoplanetary disks — stars: variables: T Tauri 



1. INTRODUCTION 

Comets are believed to be some of the most primi- 
tive bodies in the Solar System. As such, they preserve 
information about Solar nebula conditions during the 
early stages of planet formation. Samples collected from 
the comet 81P/Wild 2 during the course of the Star- 
dust mission show that the nonvolatile material of the 
comet is rich in relatively la rge grains (^few-20 /im) of 
crystalline silica tes (Brownle e et al.ll2006l : iZolenskv et al.l 
120061 : iWestphal ct al.i.2009.) . Crystalline silicates are not 
detected in the diffuse interstellar medium (ISM) and, al- 
though there is debate as to the formation mechanism of 
these grains in protoplanetary disks, they are not thought 
to form as far out in the disk as comet 81P/Wild 2 is 
believed to have originated. The Stardust results thus 
provide direct evidence for the outward transport of par- 
ticles through protoplanetary disks. Spectroscopic obser- 
vations suggest that similar processes may be widespread 
in T Tauri disks. Spitzer studies have revealed an appar- 
ent link between disk crystallinity, and grain growth and 
settling, but otherwise these disks exhibit a wide range of 
crystalline-to-amorphous silicate mass ratios and spatial 
distributions. This is suggestive of a great diversity in 
protoplanetary dust processing ()van Boekel et all 120051 : 
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iWatson et all 120091: [Olofsson et al.ll2009t) . 

Because no process has been proposed to form crys- 
talline silicates out as far as the comet-forming re- 
gions, a mechanism is needed that can explain the pres- 
ence of particles processed through high temperatures 
in cold regions of the Solar nebula. One obvious pos- 
sibility is radial mixing induced by turbulence. Tur- 
bulence is probably necessary for ga s within the disk 
to accrete (jShakura fc SunvaevI I1973D , and that same 
turbulence will result in diffusion of gaseous tracers 
and sr nall particles coupled to the gas by aerodyn amic 
forces (IMorfill fc Volkl[l98l IClarke fc Pringld [l988i ) . It 
is unclear, however, whether turbulent transport suf- 
fices to explain the observations, or whether other phys- 
ical processes are also needed. Examples of such addi- 
tional mechanisms include the ballistic launching of par- 
ticles in a wind from near the inner edge of the disk 
(jShu et al.ll2001l ). photophoretic gas-pr essure forces act- 
ing o n grains in an optica lly thin disk (jKrauss fc WurmI 
l2005t iMousis et al.l l2007f ). and radiati on pressure o n 
larger grains near the surface of the disk ('Vink ovi32009D . 

The radial transport of both gaseous and 

particle species within a turbulent disk has 

been studied by several a uthors (iGai 
Bockclee- Mor van et cH 120021: iKeller fc Gail 
Bossi 12004: iD uUcmond. Ap ai. fc WalchI 120061 
Alexander & Armit aed 12007 : iGarau 
120081 : IDullemond etall 120081: ICiesL 
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[Turner. Carballido. fc Sanol I2010D . Here, we use a 
particle-based approach to model the advection and 
turbulent transport of non-interacting dust grains 
within evolving protoplanetary disks. Our goals are to 
identify the conditions under which significant outward 
transport of particles can occur, bearing in mind the 
range of uncertainty in disk physics and evolution. We 
quantify the extent of radial transport as a function of 
the size of the grains, the initial compactness of the disk, 
the relative diffusivity of disk-gas tracers, and on the 
vertical profile of the gas fiow and particle distribution. 
Initially compact configurations that expand rapidly 
appear to have the greatest promise to explain the Star- 
dust results, though unambiguous predictions require 
accurate knowledge of the internal disk-fiow structure 
and of crystalline-silicate formation mechanisms. We 
find that the importance of the advection of solids 
within the gas flow means that the outward transport 
efficiency drops significantly for larger particles ('^few 
millimeters), and at later times, thereby limiting the 
extent of mixing uniformity achievable within the disk. 

In Section [5] we provide a summary of the relevant 
observations of silicates in both primitive Solar System 
bodies and in other protoplanetary disks. In Section |3] 
we describe our modeling approach, detailing our ID gas- 
disk evolution and particle transport models. In Section 
m we step through a set of examples designed to illustrate 
the basic effects of various model parameters. In Section 
[5l we present and analyze results for a baseline case of 
particle transport in two scenarios of radial gas flow, then 
consider the effects of varying the size of the particles, 
the diffusivity, and the initial compactness of the disk. 
We present our conclusions in Section [H] 

2. OBSERVATIONAL CONSTRAINTS 

Observations of dust in the diffuse ISM show that 
the vast majority of silicates there are amorphous, with 
sub-micron grai n sizes and a generally balanced Mg- 
Fe co mposition (jMolster fc Kempeill2005t iWooden et al.l 
|2007| ). Therefore, the crystalline silicates that we ob- 
serve in the Solar System and in disks around other 
stars are expected to have formed after the onset of 
the star-disk-system formation. Several mechanisms 
have been proposed: (1) The evaporation and conden- 
sation of silicate vapor, requiring gas heated to tem- 
peratures of T ^ 1250-1450 K (probably at less than 
an AU from the parent star 1 (lLoddersll2003HGaill 120041 : 






IWooden. Harker. fc BreaHevJ 120051: IBefl et all 119971) 77?) 
Annealing of amorphous silicate grains, requiring T > 
1000 K, (possibly occurring as far out as ^2-3 AU in very 
hot early disks) (|N uth fcJ ohnsonI 120061 : IWooden et"ar 
120071 iWestphal et al . 2009) : and (3) Shock-heating and 
annealing in disk spiral arms. The first two of these 
are equilibrium processes, whereas the third relates to 
transient events that rely on the disk being massive 
enough to produce spiral arms. It may, however, re- 
sult in the production of cry stalline silicates out as far 
as 10 AU (JHarker fc DescH [20021 : IScott fc KrotI [20051 : 
IWooden et al.l l2007[ f^. These mechanisms have differ- 
ent chemical signatures. Grain formation in long-lived, 

^ Note that although the stren gth of spiral shocks is ex pected to 
be greatest at many tens of AU IjClarke fc Lodatol 120091 ) . the gas 
in these regions is too tenuous and cold to anneal silicate grains. 



high-temperature regions of the disk is most likely to 
produce Mg-rich silicates due to the low oxygen fugac- 
ity expected to prevail in those regions. Fe-bearing and 
Fe-rich crystalline silicates likely require a water-rich re- 
gion of the disk in which to form, produced perhaps by 
migration and sublimation of icy bodies inter ior to the 
snow line, or in shocks in the outer disk (jWood cn 2008). 
They may also require transient heating mechanisms to 
form without evaporating, unlike Mg-silicates, which can 
crystallize below short-time-scale-evaporation tempera- 
tures (Nuth fc Johnson 2006). Pyroxene is thermody- 
namically favored over olivine when these minerals are 
formed by condensation, and, while annealing will tend 
to produce olivine from forsterite and pyroxene from en- 
statite, forsterite will convert to enstatite in long-lived 
(~ 10 ^ years to completion), high-temperature condi- 
tions (Gafl '20041 [Wooden. Harker. fc BrearleV 20051. 

Observationally, it has been known for some time 
that Oort Cloud comets contain crystalline-silicate mate- 
rial (i Hanner e t al.|[1994[: [Harker et al.|[2q02[: [Honda et al.l 
[200l ~Woodc n, Woodward, fc Harker[ [200411 These 
comets are believed to have formed primarily at dis- 
tances of ~ 5 — 10 AU from the Sun, in a region of the 
disk that may have overlapped that where in situ crys- 
talline silicate formation was possible. The high amounts 
of crystalline silicates recovered by Stardust, nonethe- 
less, came as a surprise, since 81P/Wild 2 is a short 
period comet, that most probably formed in the outer 
disk, a round t he curren t orbits of Uranus and Neptune 
(jWooden. Harker. fc Brearlcv 2005) . This is beyond any 
plausible source region for crystalline silicates. Subse- 
quently, observations of the comet 9P/Tempel 1 have 
confirmed crystalline silicates to be an imp ortant com- 
pone nt of that Jupiter-family comet as well (jLisse et al.l 
[20061 ). The compositional evidence from the Stardust 
samples points to diverse for mation environ ments for the 
high temperature materials (|Wooden|[20080 . A predomi- 
nance of Mg-rich silicate grains, together with the recov- 
ery of three calcium- aluminum (CAI-type) minerals that 
almost certainly formed by evaporation and condensa- 
tion at T > 1400-2000 K, suggest a substantial contribu- 
tion from the innermost, hottest, disk regions. However, 
some Fe-bearing and Fe-rich crystalline silicates were also 
recovered. There are also hints o f igneous and aque- 
ous alteration among some grain s poswiak et al.l I20T0I : 
[Stodolna. Jacob fc Lerouxl [20101 ) , which, if confirmed, 
would place additional constraints on the timing of the 
outward transport and incorporation of these materials 
into the SlP/Wild 2 cometesimals. 

Direct comparison between the Stardust Solar Sys- 
tem results and observations of other disks is difficult. 
Astronomical measurements are primarily sensitive to 
smaller orbital radii, and thus provide more impor- 
tant constraints on the degree and nature of processing 
than on radial transport. The observed mass ratio of 
crystalline-to-amorphous silicates varies g reatly, both be- 
tween disks (l^ n Bockcl ct al. 2005; Wa tson et al.|[2009[: 
[Olofsson et al.' 200 9) and likely radially within a sin- 
gle disk (l^an Bockel et al J [20041: [Olofsson et al.[ [2009). 
A study bv [Watson et al.l ()2009t ) measured crystalline 
mass fractions in the inner disks (<10 AU) of 84 clas- 
sical T Tauri stars in the Taurus-Auriga star-forming 
region. The measured mass fractions ranged from less 
than 0.5% to more than 80% despite the fact that the 
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systems were all of similar ages, 1-2 Myr, and all were 
observed within a single star-forming cluster. A signifi- 
cant correlation is observed between the crystalline mass 
fraction and the extent of dust settling toward the disk 
midplane. This may be related to the link other studies 
have found between disk crystallinity and g rain growth 
(|van Boekel et al.ll2005HOlofsson et al.ll2009li . with char- 
acteristic crystalline grain sizes of a few microns. No cor- 
relations were found relating to the mass or luminosity 
of the star, the mass of the disk, or the mass ratio of the 
star-disk system, all properties expected to affect heating 
and thermal processing within the disks. These results 
agree with a previous s tudy of Herbig Ae-star disks by 
Ivan Boekel et al.] (|2005D who found a correlation between 
disk crystallinity and stellar mass/luminosity that disap- 
peared for Mi, < 2.5Mq. 

In terms of composition, silicates within other pro- 
toplanetary disks show the same pred ominance of Mg- 
rich grains as in the Stardust samples (Ivan Boekel et all 
2004L 120051: fMolsterfc KemDedl200l I Watson et al.ll200a 



Olofsson et al.ll2009f) . There are also marginally signif- 
icant correlations between crystallinity and stellar ac- 
cretion rate. Watson et al. (2009) found that for py- 
roxene, the trend in crystalline mass fraction was in- 
verse to the mass-accretion rate onto the star, whereas 
for olivine, the trend was proportional to the accre- 
tion rate. A study bv lOlofsson et ahl (|2009[ ) found that 
disks of higher crystallinity tend to be dominated by en- 
statite (pyroxene-type) grains, and of lower crystallinity 
by forsterite (olivine-type) grains. lOlofsson et al. (20091 ) 
also suggest very heterogeneous mixing of silicate parti- 
cles in disks, reporting a higher rate of crystalline-feature 
detection for the cold (<10 AU) than the warm (<1 AU) 
spectral features. We will discuss later how some of these 
trends may be interpreted within the context of a turbu- 
lent transport model for grains within the disk. 

3. METHODS 

The simplest model for studying the radial redistri- 
bution of particles within a disk assumes that the disk 
surface density is static and that the particles are small 
enough as to be perfectly coupled to the gas. Even in 
this limit, a range of other physical effects can be im- 
portant and have been explored: mixing by (marginal) 
gravitational instability (.Boss, |2004 , 2008), the the rmal 
evolu t ion and production of anneale d silicate grains (jGaili 
120011: iBockelee-Morvan et al.) I2002D . and vertical mix- 
ing and sett ling in an MHD model that inclu des dead- 
zone effects (jTurner. Carballido. fc Sanoll2010f ). Many of 
these studies agree that turbulent mixing can, in prin- 
ciple, be strong enough to explain the presence of crys- 
talline silicates in comets. 

Our goal in this work is to systematically incorporate a 
range of additional physical processes into what is still a 
relatively simple one-dimensional model for the evolution 
of particles within a gas disk. We focus on three effects: 

(1) Imperfect coupling of particles to the gas. The 
largest particles captured by Stardust (about 
20 /Ltm) may - depending upon the gas density 
- be only marginally well-coupled. Indeed, prior 
work that includes grain size and settling effects has 
found that large grains may settle out the the mid- 
plane and experience outward advection in a 2D 



stratified disk model ()Cieslall20Q7ll2009D . Here we 
consider imperfect coupling to examine the feasi- 
bility of transporting large grain sizes radially out- 
ward into the more tenuous outer disk and comet- 
forming regions. 

(2) Disk evolution. The formation of material pro- 
cessed at high temperature almost certainly com- 
menced early in the disk lifetime, when the disk 
would have been more massive, hotter, and more 
compact than the typical T Tauri disk. As we 
will show later, the evolution of such disks can 
have a substantial impact on the radial trans- 
port of pa rticles. This has been demonstrate d ex- 
plicitlv bv iDullemond. Apai. fc WalchI ()2006D and 
iDuUemond et al.l ()2008[ ). who found that disks that 
form in initially more compact configurations pro- 
duce greater outward mixing of hot material. 



(3) 



Uncertainties in the radial flow of gas at the mid- 
plane. Although the vertically integrated flow of 
gas in an active disk is assuredly inward at small 
orbital radii, the actual magnitude (and even di- 
rection) of flow at the midplane depends upon un- 
known aspects of the angular momentum transport 
within the disk. Indeed a flared viscous disk model 
with little to n o vertical mixing yields an outward 
midplane flow (IUrpinl[l9 8l iTakeuchi fc Linl [200l 

which could be very im- 



iTscharnuter et a~ 



mmm^ 

mm, 



portant for the transpor t of particles large enough 
to have partially settled ()Keller fc Gailll2004HCieslal 
|2007, 



We describe below our implementation of these effects 
within an otherwise standard one-dimensional model for 
the evolution of a gaseous protoplanetary disk, which is 
chosen so as to be consistent with observational estimates 
of the lifetime and accretion rate evolution of protoplan- 
etary disks around low mass stars. 

3.1. Disk Profile and Evolution 

We use the thin-disk model f or the viscous evolution 
of the gas-disk surface density (!Pringlelll98lD . We also 
include loss of disk gas due to photoevaporation by the 
central star, so that the surface-density profile evolves 
according to 



5Eg 
dt 



3__d_ 
R'dR 



(R,t) 



— Sg^wind {R,t) , (1) 

where Sg is the surface density, t is time, R is radial dis- 
tance from the central star, and v is the local disk viscos- 
ity. The Sg. wind-term is the loss d ue to photo evaporation 
by E UV (jHollenbach et al.lll994D or X-ravs (lOwen et al.l 
l2Q09f ). Here we use the lAlexander fc Armita,gel (I2007D pa- 
rameterization of the results of iFont et al.l ( 20041) . which 
calculates photoevaporative mass loss due to a lO**^ s~^- 
photon fiux of EUV from the central star. We do not here 
study the evolution of particles during photoevaporation 
(long after the presumed epoch when crystalline material 
must have been incorporated into comets or their progen- 
itor bodies) ; we include photoevaporation only because it 
is essential for providing reasonable model-disk lifetimes. 
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We evolve the disk numericany using a standard 
time-exphcit-finite-differe n cing s cheme similar to that of 
iPringle. Verbunt fc Wadi (|1986D . We use 600 grid cells 
spaced logarithmically. Innermost and outermost grid 
points are placed at 0.1 and 15,000 AU, respectively. 

We assume that the disk temperature and viscosity 
are both fixed and vertically uniform. We use a tem- 
perature profile of T (R) = ToR~^^^, normalized to a 
temperature of 280 K at 1 AU. This is appropriate for 
a passive, flared disk after most of the gas infall has oc- 
curred and the initial rapid phase of accretion has slowed. 
Note that very-e arly disk temperatures are expected to 
be rauch higher (iKenvqn fc HartmannI 119871 : iBell et aD 
llQQTt iTscharnuter et al.ll2009() . However, a more precise 
treatment of disk temperature would require models or 
assumptions involving, for example, disk chemistry and 
disk-forming infall ra tes. For the disk viscosity , we use 
the a-prescription of IShakura fc SunvaevI ()1973l ). which 
gives 

ly (R) = aCsH^ , (2) 

where a is a non-dimensional scaling parameter, Cg is the 
local sound speed, and Hg is the local-disk scale height. 
Having assumed T ex _R^/^, Equation © leads to v oz R. 
We adopt a = 10~^. This value falls within the range of 
values estimated obser vationally (Hartmann ct al. 1998; 
IHueso fc Guinotl[20()5l : lKing. Pringle. fc Livioll2007i 'l and 
allows dissipation of our model disks in less than 10 Myr. 
The parameterization in Equation ([2]) assumes that the 
source of the viscosity is turbulence within the gas disk. 
Equation ([2]) is also important for establishing the tur- 
bulent diffusion properties for the particle ensemble, as 
discussed in Section [331 

The initial surface density distribution of the gas disk 
depends on the mass and angular momentum distribu- 
tion of the p arent cloud that collapsed to form the star - 
disk system (|Lin fc Pringlel[T990HHueso fc Guillotll2005l) . 
Our evolving-disk models do not use a steady-disk pro- 
file, but rather assume a finite disk, so that at t = 0, the 
surface density is. 



Eg(i?,t = 0) = 



Mo 
Sniy 



1- 




cxp 



R_ 

Rd 



(3) 



where Mq is the t = accretion rate onto the central 
star, i?in is the inner-disk boundary (set equal to the 
inner-grid boundary, R1/2 = 0.099 AU), and i?d is a 
variable controlling the compactness of the t = profile. 
In our nominal-disk model, Rd — 20 AU, and the ini- 
tial disk mass of 0.03 Mq dictates Mq = 1.8 x 1O~''M0 
yr~^. The evolution of this model is shown in Figure [TJ 
As the disk evolves, the surface density drops, and the 
disk spreads outward before dissipation via photoevap- 
oration occurs after ~5.5 Myr. That disks may form 
in initially compa ct co nfigurations has been discussed 
by several authors (iLin fc Pringle 1990; Hucso & Guillot 
120051 : iDullemond. Natta. fc Testil 1200 61 and may be in- 
directly supported by the ste e p disk surface densities 
calculated by Wcidens chillind (|1977b( ) for estimates of 
the minimum -mass Solar nebula, and, more recently, by 
iDeschI ()2007[ ). Depending on the strength of the disk vis- 
cosity, a compact initial configuration may also be nec- 
essary to process expected disk masses within observed 
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Fig. 1. — Surface-density profiles for the disk gas. The purple 
curves track the time evolution of the nominal-disk model. The 
yellow curve is the t = profile for the most compact (R^ = 5 
AU) disk model. The thick grey curve is the profile of the static 
steady-disk model. 

disk lifetimes. 

In Section WM we vary R^ (retaining Md.o — 0.03AfQ) 
to consider the effects of varied disk-forming conditions 
on the transport of particles. To facilitate a compari- 
son with other models, we also consider two models of 
dust motion within a static disk, including a steady-disk 
model where i?d = 00 (retaining Mq — 1.8 x lO^^Af© 
yr~^). As shown in Figure [l] this steady-disk model 
results in a disk profile that is matched to our i = 
nominal-disk profile in the inner-disk regions, but main- 
tains a shallower profile and much larger gas surface den- 
sities in the outer disk. 

3.2. Two Cases for Disk Gas Radial Velocity 

The mean radial flow of gas experienced by solid par- 
ticles within the disk depends on aspects of disk turbu- 
lence that cannot yet be predicted from flrst principles. 
We maintain this uncertainty explicitly by considering 
two bounding cases for the gas flow, the "accretion-flow" 
case and the "midplane-flow" case. For the accretion- 
flow case, we use the vertically averaged gas velocity, 
^acc, (|Pringlel ll98in 



d 



Sgi?i/2 ai? V " s 



vll^R^/^ 



(4) 



This velocity can be derived from the ID fluid equa- 
tions and results in a gas flow that is predominantly in- 
ward. Importantly, however, there is a region of outward- 
flowing gas in the outer disk where the disk is expanding. 
The boundary of this region moves outward as the disk 
evolves and spreads. 

In a real two-dimensional (2D) disk, the radial flow 
may vary with height above the disk midplane, and par- 
ticles that settle toward the midplane may experience a 
very different radial-gas velocity than that predicted by 
Equation Q. In particular, 2D {R,z) disks in which 
the viscosity is given by a generalized a model with 
no vertical mixing of angular momentum yield a n out- 
ward flow at the midplane at most disk radii (lUrpinl 
119841: iRozvczka. Bodenheimer fc Belli Il994[ ). This flow 
is highly conducive to the transport of inner-disk grains 
to the outer disk For the midplane-flow case, we follow 
iTakeuchi fc LinI ()2002[ ) and use the azimuthally symmet- 
ric Navier- Stokes equation to calculate Umerid, the gas 
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Fig. 2. — Radial advection velocities of the disk gas (black/grey) 
and 20 fim dust (red/pink). The black and red curves are for 
the accretion-flow case, and the paler-colored curves are for the 
midplane-flow case. Solid lines plot t = 0, and dashed lines plot 
t = 1 Myr, both for the nominal-disk model. Dotted lines plot the 
gas velocity for the steady-disk model. 

velocity at the midplane: 
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dz 



R' 



,drtg 

~dt 



(5) 



Here, pg is the local-gas volume density and z is height 
above the midplane. ilg is the local orbital angular veloc- 
ity of the gas, determined via fo rce balance in the radial 
direction (JTakeuchi fc Linll2002| ). 



mi 



GM^R 

(i?2 + ^2)3/2 



J_5p 
PgdR 



(6) 



In this equation, M-^ is the mass of the central star (equal 
to 1 Mp,), and p i s the local gas pressure. Note that while 
iTakeuchi fc LinI ()2002|) assumed a power law for the gas 
distribution and used that to derive a simplified expres- 
sion for the flow structure of the gas, we have used the 
input of the disk surface density from our disk-evolution 
model to solve Equation ([5]) numerically by assuming 
vertically uniform temperature and viscosity. 

Figure [2] shows the two gas- velocity cases at t = and 
at i = 1 Myr. Note that the outward expansion of the 
disk dominates both velocity cases at large disk radii. 

3.3. Particle transport 

Our model for particle transport tracks an ensemble of 
particles within a ID gas-disk environment and subjects 
the particles to two non-gravitational forces: (1) aerody- 
namic drag against the mean gas flow and (2) turbulent 
diffusion of the ensemble within the disk gas. Coagu- 
lation is neglected. Our particle trajectories therefore 
consist of two components: (1) mean radial motion as- 
suming Epstein drag against the mean-disk flow and (2) 
a random-walk component due to turbulent diffusion. 

These two physical effects are represented numerically 
by two components to the dust radial velocity: Vsrd, the 
mean radial velocity induced by Epstein drag, and Wturb, 
the turbulent random walk velocity. These velocities are 
calculated at each evolutionary time step of the gas disk 
model for each grid point of the model disk; values be- 
tween grid points are linearly interpolated as necessary. 



Together, these velocities are used to integrate the dust- 
grain trajectory as: 

Td (t + At) = ra (t) +AtX {Vsrd ± "turb) , (7) 

where r^ is the dust-particle radial position in the disk, 
and At is the time step. The time step is chosen locally, 
to insure accurate integration of particle trajectories, but 
is capped at the global disk evolution time step. The fi- 
delity of this method in reproducing analytical results 
for trace-particle diffusion and transport in a static disk 
is discussed in the Appendix. The relative motions of 
the dust and gas can vary from tightly coupled to diver- 
gent depending on the surface-area-to-mass ratio of the 
simulation particles and the local density of the disk gas. 

3.4. Advection 

We solve for the mean radial velocity of the dust grains, 
Vsrd, subject to Epstein drag against the mean gas flow. 
In the Epstein-drag regime, the size of the particle is less 
than the mean free path of a gas molecule, and the drag 
scales with both the particle velocity and t he thermal 
velocity of the gas as (|Weidenschillingin-977a|) : 






Vd 



(8) 



where Cr and m^ are the surface-area-to-mass ratio and 
mass of the dust grain, respectively, pg is the local gas 

density, Wthorm — \/o/ttCs is the local thermal velocity of 
the gas, and Vg and v^ are the gas and dust velocities. 
In our model, we assume that Pg = Pg.mid = '^g/V^nHg, 
the gas density at the disk midplane. In the accretion- 
flow case, inputting the midplane density to Equation ([8]) 
is an approximation that selects for the tightest possible 
coupling between the gas and the particle ensemble at 
a given R; it leads to the most conservative estimate 
of outward mixing in this gas-flow case. However, the 
midplane density is the most appropriate choice for use 
in the midplane-flow case, which presupposes that the 
entire particle ensemble has settled to the disk midplane. 
To solve for the steady mean radial velocity of the dust 
grains, we must include both the drag force and the or- 
bital motions in a polar coordinate system. (See the ap- 
pendix for the expli cit force- balance e quati ons of particle 
motion.) We follow iTakeuchi fc LinI (|2002l ) and simplify 
these forces by assuming that the radial acceleration of 
the grain is zero, dvr^d/dt w 0, and that the azimuthal 
acceleration corresponds to a change in the Keplerian 
velocity with R, dv^fi^d/dt sa — UKWr,d/2i?. The equations 
for the radial and azimuthal motions of a grain now pro- 
duce two equations for the steady mean radial velocity, 
Vsrd, as a function of its steady azimuthal velocity, Vs^d- 



'^ srd '^r,g ~r 
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(9) 



which we solve iteratively for Vsrd using the gas-disk bulk 
flow from our disk model at each radial grid point. In 
the Appendix, we demonstrate that the trajectories pro- 
duced using the gridded-Ws^d values agree well with di- 
rect numeric integrations of the force-balance equations 
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of grain motion for particles in disks with Epstein-drag 
stopping times ranging from less than 10^^ Kepler times 
to greater than ten Kepler times. 

Figure [2] shows how the mean radial velocities of the 
dust compare to the radial velocity of the gas. Very small 
grains are well coupled to the gas flow, but the crystalline 
silicates larger than 10 /im associated with the Stardust 
collection are large enough that at large distances from 
the star (tens of AU) and at late times in the disk evolu- 
tion, the dust motion significantly departs from the gas 
flow. 

3.5. Turbulent Diffusion 

To model the turbulent diffusion of the particle ensem- 
ble, we add a random walk to the individual grain mo- 
tion via the addition of a turbulent velocity component, 
^turb- This turbulent velocity is set by the dust particle 
diffusivity, D^ , and by the time step of the random walk 
motion; it is given by 



I'turb = ± 



2i?p 

At 



(10) 



In our particle-transport model, the time step of the ran- 
dom walk is limited by the time step of the gas-disk- 
evolution model, but is also allowed to be no greater 
than the local Kepler time, l/fix (^k is the Keplerian 
angular velocity), and no greater than either ARc,/vsrd 
or Ai?5/wturb, where AR^ is the radial distance across 
five grid cel ls in the direction of partic le motion. 

We follow I Youdin fc Lithwickl ()2007[ ) in calculating the 
radial diffusion coefficient of the dust grains via 



Dp^D, 



(l + ri) 



(11) 



where Dg is the gas diffusion coefficient, and Ts = 
^K^top- Here, istop is the exponential-stopping time. 
In the Epstein-drag regime, it is given by igtop = 
3/ (CiiPgUthcrm)- The dust and gas diffusion coefficients 
are equal throughout most of the disk, except in the out- 
ermost regions where the gas surface density is very low, 
and the dust diffusion coefficient drops to zero. In our 
fiducial models, we assume that because the disk viscos- 
ity is derived from turbulence within the disk, Dg — v. 
In section [Ol we vary Sc = vjD^ to evaluate the im- 
pact of this assumption on the potential for the outward 
mixing of inner-disk grains. 

While Equation ([T0| provides an accurate represen- 
tation of the particle diffusivity, it is not sufficient for 
producing accurate diffusion of the particle ensemble in 
our model gas disk. Among other things, it does not 
account for the mass distribution of the disk gas. In a 
real random walk, the time and distance of each step 
vary according to the local properties of the gas envi- 
ronment. In our model, however, the time step is fixed 
by other, mostly numerical, considerations. Therefore, 
to allow our particle ensemble to diffuse according to the 
actual gas distribution, we must calculate an appropriate 
probability for whether a given particle in the ensemble 
will step radially inward (pin) or outward (pout)- 

The proper weighting between pi„ and "Pout is a func- 
tion of the imposed time step. Instantaneously {At — )• 0), 
we know that Pin/Pout = 1, because At = repre- 
sents only a single random walk encounter. However, 



for At > 0, the probability that multiple encounters will 
occur within that time step becomes nonzero. Further- 
more, additional encounters will occur not at the ini- 
tial position being considered, but either inward or out- 
ward of that location depending on the outcome of past 
encounters. Therefore, the properties of the diffusive 
medium (the gas) outside the local point of interest play 
a role in determining Pm/Pout for At > 0. The spatial 
extent of this region of interest is set by how far particles 
may travel for a given At. For the particle ensemble, 
this is the root-mean-square of the displacement of the 
diffusing particles, or 



AR 



(Ax)' 



vtu^At ^ y/2DpAt . (12) 



The weighting for pin is set by the gas properties be- 
tween R — AR and R; and for pout by the properties 
between R and R + AR. Asymptotically, the diffusion 
must yield everywhere a uniform particles(dust)-to-gas 
ratio. Therefore, after an infinite time, the region with 
more gas mass must also have proportionally more dust 
particles. However, because regions of higher diffusiv- 
ity reach the steady state of uniform concentration more 
quickly, for a finite time step, proportionally more parti- 
cles will also mix into a region of higher diffusivity than 
into one of lower diffusivity. Therefore, 



Pin _ (MgXDp)^ 

Pout (Mg X Dp)^, 
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DpT^gRdR 
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(13) 
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To evaluate Equation (J13p within our model gas disk, 
we assume uniform, mean values of the gas surface den- 
sity and particle diffusivity between each grid point. Be- 
cause we calculate pi„ once at each disk-evolution time 
step, but At of the particle motions may vary locally 
for the individual grains, we calculate pin at each grid 
point for both Aicvoivo and ^Atcvo\■vc^ then interpolate 
parabolicly for other time-step sizes as needed, using 
Pin {At = 0) = 0.5. 

We have verified that this random-walk method does 
a good job of reproducing the expected diffusion pro- 
file of a contaminant within a ID gas disk. Figure [3] 
presents a comparison of the distribution of particles 
using our particle-transport model compared to the ex- 
pected distribution for an analytical example derived by 
[Clarke fc Pringld (|1988[ ). A more detailed discussion of 
the fidelity of our particle-diffusion model may be found 
in the Appendix. 

4. EFFECTS OF INDIVIDUAL TRANSPORT PROCESSES 

To illustrate how each of the model components in- 
fluences particle transport, we first consider a set of re- 
duced models of successively greater complexity. We be- 
gin with a diffusion-only case for perfectly coupled par- 
ticles, which we run keeping the gas disk profile fixed at 
its nominal t — {) form. We then add in the effects of 
radial gas velocity, particle size, and disk evolution. Un- 
less otherwise stated, these simulations use 1000 particles 
initiated at 1 AU in the nominal-disk model. 

Figure |4] displays an image of the relative number of 
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Fig. 3. — The fraction of particles per radial bin for our particle- 
transport simulations (pale lines) versus the analytical-expected 
values (dark lines). This figure corresponds to th cClarke fc Pringld 
II1988I ) example where u oc S^ oc R^ in a static disk with an 
accretion- flow velocity, tiacc = —l.bv/R, and y = 4.941 X 10^* 
cm'^ s~^ at 1 AU. The simulation was initiated at t = with 10* 
particles at 40 AU. 
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Fig. 4. — Basic diffusion: Map of the number of particles per 
grid space as a function of radial position and time. We initiated 
1000 well-coupled particles at 1 AU (denoted by the green line) 
in a static disk with zero radial-gas velocity. Note that the radial 
width of the grid cells increases exponentially with distance from 
the host star. 

particles in each radial bin over time for the diffusion- 
only siniulatior{j. Diffusion results in rapid (within ~ 
10^ yr) radial spreading over tens of AU, such that the 
particle distribution closely approximates that of the gas 
within a fraction of the disk lifetime. Some particles are 
also lost past the inner edge of the simulated disk. 

Next, we include the effects of gas velocity by using 
the two radial-gas- velocity cases outlined in Section 13.21 
The distribution maps for these simulations are shown 
in Figure [5] (al, a2), and here the gross effects of mostly 
inward vs mostly outward gas velocities are clear. In 
the accretion-flow case, where most of the gas flow is 
inward, the majority of the particles are quickly lost onto 
the central star. Only ^2% of the particles remain in 
the simulation after 10'* years. However, the particles 
that do remain in the disk reach the outer-disk region 

* Note that the bin size increases exponentially with R, so this 
representation looks somewhat different from a plot of the particle 
surface density or concentration. 



and continue to move outward, because the steep gas 
distribution in the outer disk causes the accretion flow 
there to be outward. In the midplane-flow case, less than 
10% of the particles are lost. Particles are instead swept 
rapidly outward with the gas flow and become trapped 
at the edge of the static disk-gas distribution, where the 
value of the disk surface density falls rapidly toward zero 
(at about 340 AU). 

After turbulent diffusion and gas velocity, we add in 
the effects of particle size. The results of simulations us- 
ing 20 /im particles are shown in panels (bl) and (b2) 
of Figure [5] . The effects of particle size are most clearly 
visible in the midplane-flow case where particles are still 
advected outward as before, but pile up at smaller ra- 
dial distances (around 90 AU) than in the simulations 
with the well-coupled particles. What is happening here 
is that beyond ^115 AU, the disk gas is too diffuse to 
support the outward advection of particles of this size. 
Instead, headwind drag dominates particle motion, push- 
ing particles inward. In the accretion-flow case, the ef- 
fects of a non- negligible particle size are similar. Though 
the dominant behavior of particle motion here remains 
the loss of most of the particles inward onto the parent 
star, those 20 /im particles that do reach the region of 
outward accretion flow are also confined to smaller radii 
than their massless counterparts. 

Finally, we include the effects of disk evolution. As the 
disk evolves, the surface density drops. Panels (cl) and 
(c2) of Figure [5] show that in the accretion-flow case, 
this drop in density results in a sustained inward loss 
of the particles; in the midplane-flow case, the majority 
of particles are still retained and advected outward. In 
fact, some particles reach greater distances than in the 
static-disk model, because, as the disk expands, the sur- 
face density no longer drops off as steeply around 100 
AU. However, after about t = 0.5 Myr, the decrease in 
disk surface density begins to dominate over the effect 
of the disk expansion. The gas can no longer support 
20 /xm particles at such large distances, and the particle 
distribution begins to move inward. 

So far, we have presented simulations of particles ini- 
tiated quite close to the parent star (i? = 1 AU). In 
panels (dl) and (d2) of Figure [S] we see that if the par- 
ticles are initiated further out in the disk, the particle 
motions are roughly the same, with the major excep- 
tion that significantly more particles are retained in the 
accretion-flow case. For particles initiated at 10 AU in 
this case, about 60% remain in the disk after 10^ years. 
In the midplane-flow case, particles initiated at 10 AU 
reach the furthest regions of the disk slightly faster, but 
otherwise behave the same as when initiated at 1 AU. 
Therefore, while both gas-flow cases can turbulently dif- 
fuse particles over wide regions of the disk, the accretion- 
flow case selectively retains particles initiated at larger 
radii and rapidly loses particles initiated at smaller radii. 

To summarize, the basics of particle mixing in these 
transport simulations are that: 

1 . Turbulent diffusion rapidly spreads the particle dis- 
tribution over many tens of AU. 

2. Where advection tends toward the loss of particles 
inward onto the parent star, there will be a pref- 
erential rapid depletion in grains initiated close to 
the parent star. 
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Fig. 5. — Diffusion scenarios: Maps of the diffusion of particles as in Figure 4 where the green Une denotes starting location. The right-side 
panels all include gas velocity in the accretion-flow case, and the left-side are for the niidplane-flow case. The simulations in (bl) and (b2) 
add in the effects of non-zero particle size. In (cl) and (c2) the disk surface density is also evolved with time, and the (dl), (d2) simulations 
are just like those above, but with the particles initiated at 10 AU. 
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3. Even in the case of outward-flowing disk gas (the 
midplane-flow case), small dust grains may still spi- 
ral inward at radii (or epochs) where the disk is 
locally tenuous. 

All of these effects will be reflected in the availability of 
inner-disk particles in the comet-forming regions. 

5. RESULTS 

In this section, we study how the parameters of the 
system affect the efficiency of mixing inner-disk parti- 
cles out to large distances. These parameters include the 
structure/direction of the radial gas flow, the sizes of the 
dust particles, the diffusivity relative to the disk viscos- 
ity, and the initial distribution and evolution of the disk 
gas mass. We vary these parameters relative to a fidu- 
cial model that assumes 20 /im particles (all particle sizes 
assume a dust-particle internal density of 3 g cm~^), a 
Schmidt number of S'c = 1 {Dg ~ v\ and an evolving 
nominal-disk model. The nominal-disk begins with a to- 
tal gas mass of 0.03 Mq and an exponential fall-off radius 
of Rd = 20 AU, so that at t = 0, 61% of the disk mass is 
concentrated within 20 AU of the star. 

In each model, the particles (10,000 for the accretion- 
flow runs, 2,000 for the midplane-flow) are initialized 
with a uniform radial distribution between 0.5 and 10 
AU. Although this is not an exact match to the ini- 
tial gas distribution at these radii, it allows us to fully 
sample the likely source regions for crystalline silicates 
within the disk. We are interested, primarily, in the 
question of what fraction of these particles reach the 
"comet-forming region" (deflned, for this paper, as the 
disk beyond 25 AU) at any time during the evolution of 
the disk. However, in some physical models, crystalline 
silicate formation is confined to smaller radial extents. 
Therefore, we also consider a subset of our simulation 
particles divided into three smaller source regions: the 
inner-quarter region (0.5-2.5 AU, 2106 particles in the 
accretion-flow simulations, 421 particles in the midplane- 
flow simulations), the inner-half region (0.5-5 AU, 4737 
and 947 particles), and the outer-half region (5-10 AU, 
5263 and 1053 particles). 
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Fig. 6. — Fraction of simulation particles beyond 25 AU for dif- 
ferent source regions compared to the fraction of the total disk 
gas mass currently beyond 25 AU (divided by three). Pale-colored 
lines denote particles well-coupled to the gas; dark colors denote 
20 /im-sized particles. The inner-quarter source region is 0.5—2.5 
AU; the inner-half source region is 0.5-5 AU; the outer-half source 
region is 5-10 AU. Simulation run in the nominal disk model using 
the accretion-flow case. 
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Fig. 7. — No rma lized concentration of particles beyond 25 AU. 
[Sec Equation I I14I I.] The pale-colored lines denote well-coupled 
particles. The dark lines are for 20 /xm particles. Nominal disk 
model. Accretion-flow case. 



5.1. Transport in Different Radial Gas-Flow Cases 

Figure [6] shows the baseline results for particle mixing 
if the effective gas velocity seen by the particles is that 
of the accretion-flow case. We plot the fraction of 20 /im 
particles beyond 25 AU that originate from each of the 
source regions, together with the same quantity for per- 
fectly coupled particles (i.e. of negligible size). We also 
show the gas fraction beyond 25 AU. As expected from 
the reduced models, diffusion is fast enough to mix some 
particles upstream into the comet-forming region. In- 
deed, the first particles pass 25 AU in less than 20,000 
years. However, loss of particles by accretion onto the 
star is also rapid, especially close to the inner edge of the 
disk. After 10^ years only ^30% of all 20 /im particles 
remain in the disk, and most of those originated in the 
outer half of our total source distribution (the peak frac- 
tions of particles in the comet-forming region are 2.6% 
from the inner-half zone, but 11% from the outer-half 
zone). The lifetime of particles that do manage to at- 
tain large radii is also limited, first by the fact that the 



transition radius that divides gas inflow from outflow it- 
self moves out as the disk evolves, and second because 
declining gas densities in the outer disk eventually pre- 
clude retention of 20 fim particles. These effects mean 
that by i = 10^ years, less than 1% of 20 /im particles 
exist in the disk beyond 25 AU, even if they started in 
the outer-half source region. By i = 2.1 Myr, all have 
been lost inward onto the parent star. 

To quantify the effect of these mixing patterns on the 
abundance of hot, inner-disk grains in the comet-forming 
regions, we plot in Figure [7] the normalized concentra- 
tion, Cn, of source-region particles beyond 25 AU as a 
function of time. We deflne the normalized concentra- 
tion as the number of particles in a region divided by the 
gas mass in that region, then normalized by the initial 
number of source-region particles divided by the t = 
gas mass in that source region: 
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Fig. 8. — Midplane-flow case: Normalized concentration of sim- 
ulation particles beyond 25 AU for all particles initiated between 
0.5 and 10 AU. Nominal disk model. 

With this definition, Cn is hke a sealed dust-to-gas ra- 
tio, where Cn — 1 means that the ratio of particles-of- 
interest-to-gas in the region of interest (the outer disk) 
is the same as where those particles originated at the 
time that they formed. We note that while the denom- 
inator of Equation (ITi|) is almost constant across our 
source region (i.e. the initial dust to gas ratio is a fixed 
value to within 10%), the normalized concentration be- 
yond 25 AU is an average over a broad region. The local 
values can vary widely, depending upon the parameters 
of the disk. For the accretion-flow case simulations, the 
highest local values of the normalized concentration typ- 
ically lie just outside of 25 AU, with a lower local Cn in 
the regions hundreds of AU from the parent star. 

From Figure [71 we see that the peak values of the nor- 
malized concentration - and thus the maximum extent 
to which the outer disk can be contaminated by particles 
from the hot inner regions - is relatively modest. For 
20 /im particles, the peak in Cn from the all source re- 
gion is approximately 6%. This is actually higher than 
the value obtained from the outer half source region, re- 
flecting the fact that although a greater fraction of parti- 
cles mix outward from this region, there is less mass there 
to start with. We also observe that particle size strongly 
influences the lifetime of particles in the outer disk. For 
well-coupled particles, Cn declines only modestly from 
its peak at a few hundred thousand years out to several 
Myr. This reflects the fact that well-coupled particles 
that survive the initial evolution have thoroughly mixed 
with most of the entire disk gas mass, so that their sub- 
sequent loss is largely in proportion with the disk-gas ac- 
cretion rate. By contrast, 20 fim particles are lost much 
faster, and only signiflcantly contaminate the outer disk 
between '-^ 5 x lO"' and 8 x 10^ years. 

Dramatically higher efficiencies of particle retention 
are obtained if, instead of experiencing the mean gas flow, 
the particles settle and experience an outward fiow at the 
disk midplane (recall that this requires little or no verti- 
cal mixing of angular momentum, so that the steep radial 
density gradient at the midplane results in outflow) . Re- 
sults for this midplane-flow case are plotted in Figure [51 
In this limit, few 20 //m particles are initially lost inward 
onto the parent star. 98% remain in the disk after 10^ 
years, and around 90% of particles from each source re- 



gion are transported to beyond 25 AU, with only a slight 
selection against particles from the innermost regions of 
the disk. This means that the normalized concentration 
curves for the different source regions are all nearly the 
same, only scaled relative to a given source region's in- 
nate ability to contaminate the outer disk - the relative 
fraction oi t — disk mass in that source region. 

Because 20 fim and well-coupled particles are swept to 
large distances and remain there for a long time, there 
is a period when the change of normalized concentra- 
tion beyond 25 AU is a function solely of the loss of disk 
gas mass via accretion. Without the concurrent loss of 
accreting particles, the particle-to-gas ratio reaches val- 
ues even higher than in the original 0.5-10 AU source 
regions. This is particularly so for the well-coupled par- 
ticles, which never leave the outer disk once they have 
entered it. The normalized concentration of 20 fim par- 
ticles drops to half-maximum at around t ~ 2.6 Myr, 
about three times the contamination time-window given 
in the accretion-flow simulations. 

Note, however, that our model cannot account for the 
radial mixing that would occur in the presence of verti- 
cal mixing and a vertically stratifled gas flow. If the gas 
flow is rapidly inward far from the midplane and outward 
near the midplane, then vertical mixing will lead to ra- 
dial raixing, even in the absence of radial diffusion. iCieslal 
(|2007l [20091) discusses in detail how a given degree of ver- 
tical mixing can lead to enhanced mixing of grain popula- 
tions or to the partial segregation of inward-flowing and 
outward flowing populations. Therefore, the normalized 
concentration for the midplane-flow case merely repre- 
sents the extreme upper limit on the extent and dura- 
tion of the contamination of the outer disk by inner-disk 
material. 

The results imply that the efficacy of the fiducial disk 
model (Rd — 20 AU, Sc = 1) for outward transport 
of 20 fim grains depends upon both the location (and 
extent) of the the source region for crystalline silicates, 
and upon the radial gas flow experienced by the parti- 
cles. In the accretion flow case, the peak of contami- 
nation of the outer disk (beyond 25 AU) would result 
in a crystalline-to-amorphous ratio of ^-^6%. Attaining 
this limit requires that at i = all silicate grains in the 
entire inner disk out to 10 AU were crystalline (and like- 
wise assumes that all beyond that were amorphous). If, 
instead, the source of crystalline silicates is conflned to 
inward of 2.5 AU then the peak Cn beyond 25 AU is only 
about 0.2%. Moreover, we flnd that even these modest 
degrees of contamination are relatively short-lived. 20 
fim-sized grains (matched to the upper-size end of the 
Stardust grains) only have a contamination lifetime in the 
outer disk of < 1 Myr. Unless the timing of the epoch 
when grains became assembled into cometary material 
happened to coincide with the peak of the contamina- 
tion, the crystalline fraction would be diluted, either by 
grains grown locally in the cold comet-forming region or 
by those raining inward from the even-colder outer re- 
gions of the disk. Finally, our results for the accretion 
flow case might be yet further reduced by the fact that 
some crystalline silicate formation mechanisms may not 
yield 100% crystalline fractions out to large radii. For 
example, if the primary source of crystalline silicates is 
shocks out to ~10 AU, but a sizable fraction of silicates 
in that region remain amorphous, then that cuts into the 
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Fig. 9. — The mean-radial velocity of particles of different sizes 
in the nominal disk model at i = 0, accretion-flow case. Sizes listed 
assume an internal-particle density of 3 g cm~^. 

potential maximum in the crystalline-to-amorphous ratio 
produced in the outer disk. 

In the midplane-flow case, on the other hand, substan- 
tial outward transport can occur almost irrespective of 
the location of the source region. If there is an outward 
flow of gas at the disk midplane, then grains that are 
sufficiently w ell-settled will readily reach large distances 
(|Cieslall2009fl . An outward-midplane flow also allows ma- 
terial that has been transported outward to persist there 
for a longer period of time. Our midplane-flow simula- 
tions produce an upper-limit maximum Cn in the outer 
disk of 30% for the inner-quarter source region. Even if 
vertical mixing means that only a fraction of the maxi- 
mum can be achieved, it is still plausible to believe that 
the CAI-like grains captured during the Stardust mis- 
sion may have reached the outer disk via mixing within 
the disk gas alone. This sc enario may also be c onsis- 
tent with the observations of I Watson et al.l (12009') , who 
flnd higher crystalline-to-amorphous silicate mass ratios 
in disks where the grai ns appear more s ettled toward the 
disk midplane, and of lOlofsson et al.l ()2009[ ). who flnd 
that possibly a large fraction of crystalline silicates exist 
in the outer regions (;$10 AU versus <1 AU) of some 
disks. 

5.2. Transport of Differently Sized Particles 

In Section |4l we demonstrated that non- negligible par- 
ticle size bars particles from the outermost regions of 
the disk, because at large distances, the gas is so ten- 
uous that the particles decouple from the gas motions, 
and inward-pointing head- wind drag dominates. Figure 
[3] plots the dust- mean-radial velocities for several differ- 
ent particle sizes in the nominal-disk model at i = 
for the accretion-flow case. (In the midplane-flow case, 
the curves are similar, only shifted toward outward ve- 
locities as in Figure [2]) Though turbulent diffusion can 
send particles upstream of their average radial- velocity 
flow, grains are effectively barred from the regions of the 
disk beyond where their average-radial velocity falls to 
large, inward-pointing values. It is clear from Figure 
[9] that the range within the disk that dust grains may 
occupy depends strongly on the particle size, and that 
the mixing achievable out to the comet-forming regions 
will become sharply limited as one approaches millime- 
ter particle sizes. Also, where this cutoff occurs changes, 
not only in time as the disk thins, but also from disk to 
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Fig. fO. — Maximum values of normalized concentrations beyond 
25 AU as a function of grain size. Solid lines denote simulations run 
in the accretion-flow case; dashed lines with black markers denote 
the midplane-flow case. Error-bars mark the \/N noise of the peak 
values. Nominal disk model. Grey line marks the baseline case of 
20 /^m-sizcd particles. In the midplane-flow case for the smallest 
particles, peak values shown are for t = 5.2 Myr, just prior to 
photoevaporative dissipation of the disk. 

disk depending on the total mass and mass distribution 
of each system. 

We parameterize the outward mixing of grains as a 
function of size by plotting in Figure [10] the peak values 
of the normalized concentration of grains beyond 25 AU 
as a function of grain size. The error bars show statistical 
errors based on the number of particles beyond 25 AU 
in each run. We plot the maximum values of Cn from 
each of the different source regions, as well as the upper- 
limit values provided by the midplane-flow simulations 
(dashed lines with black markers). While the magni- 
tude of outward mixing varies by orders of magnitude 
across these populations, the trends with grain size are 
the same across all sets: as seen previously, the highest 
normalized concentrations of inner-disk particles in the 
outer disk occur for the smallest grain sizes, up to a few 
tens of microns, and the peak-CAr values drop off s harply 
at around m illimeter grain sizes. The results of iCieslal 
(|2007l l2009i) for a vertically stratifled disk flow suggest 
that the loss of larger grain sizes may be partially miti- 
gated by settling toward the midplane, thereby increas- 
ing the relative outward transport and retention of these 
particles. Still, the concurrent drop in our midplane-flow 
upper-limit concentrations suggests that this can only 
hold back the loss of large grains from the outer disk to 
a point. For 2 mm-sized particles, for example, even the 
midplane-flow runs yield a peak normalized concentra- 
tion beyond 25 AU for all simulation particles of only 
-3.4%. 

The other important variable to consider is the tim- 
ing of these mixing events. In Figure [Til we depict, as 
a function of particle size, the time-frame over which 
the normalized concentration of inner-disk particles in 
the comet-forming region is at a value of half-maximum 
or greater. For each source region (for the accretion- 
flow case simulations) we plot the latest time of half- 
maximum concentration prior to the peak and the ear- 
liest time of half-maximum after the peak. (Due to 
limited particle statistics, some of these quantities are 
quite noisy.) We also include the boundary of post-peak 
half-maximum concentration in the midplane-flow sim- 
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Fig. 11. — Time- windows for the peak values of normalized con- 
centrations beyond 25 AU as a function of grain size and for dif- 
ferent source regions. The lower set of solid curves mark the latest 
pre-peak time of half-maximum concentration in the accretion-fiow 
case. The upper curves mark the earliest post-peak time of half 
maximum concentration for the accretion-flow case (solid lines) and 
the midplane-flow case (dashed lines with black markers) . Nominal 
disk model. In the midplane-flow case, the smallest particles re- 
main in the outer disk longest and never drop below half-maximum 
concentration beyond 25 AU. 
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Fig. 12. — Normalized concentration beyond 25 AU at t = 1 Myr 
as a function of grain size. The largest particles are all depleted 
from beyond 25 AU by 1 Myr. Dashed lines with black markers 
denote the midplane-flow simulations. Nominal disk model. 

ulations to illustrate when headwind-drag inspiral will 
overwhelm other physical effects even in the most ex- 
treme gas-flow scenario. 

Figure [IT] demonstrates that the smallest inner-disk 
particles remain mixed into the outer disk for millions of 
years, and that even inner-disk grains as large as 20 yum 
may be noticeably present at large distances for ~l-2 
Myr. However, while the mixing of all grains outward 
is relatively rapid in these simulations, not only do the 
larger particles experience relatively less outward mix- 
ing, but that mixing is similarly short-lived. Even in 
the midplane-flow simulations, post-peak half-maximum 
concentration is reached by i ^ 6.2 x 10^ years for 0.2 
mm particles and as soon as i ~ 1.1 x 10 years for 2 
mm particles. The fact that larger inner-disk particles 
are short-lived in the outer disk is demonstrated more 
bluntly in Figure 1121 where we plot the normalized con- 
centration as a function of grain size at t = 1 Myr. While 
the mixing effects are relatively uniform for small parti- 
cles at this time, no particles 0.2 mm and larger remain 



in the outer disk in the accretion-flow simulation, and 
the upper-limits of the midplane-flow simulations suggest 
that even in that extreme those large-grain populations 
are in rapid decline, if not already disappeared. 

These results imply that it is the upper envelope of par- 
ticle sizes recovered by Stardust that provides the most 
stringent constraints on disk evolution and particle trans- 
port. In our simulations, particles at least as large as 2 
/um appear to behave as though well-mixed with the disk 
gas. Therefore, the presence of high-temperature parti- 
cles in this size range within a Jupiter-family comet can 
be largely explained by outward mixing of disk mate- 
rial aided by some disk expansion from an initially com- 
pact state. However, the particles recovered by Stardust 
include crystalline silicates as large as 20 /im. While 
the outward transport of 20 /im particles to the comet- 
forming region does occur in our models, these parti- 
cles do not remain well coupled to the disk gas indefi- 
nitely. This suggests that either (1) mixing of particles 
and planetesimal formation in the Solar nebula occurred 
in a disk massive enough for 20 /xm particles to remain 
well coupled to the gas motions, or that (2) the forma- 
tion of comets or cometesimals occurred within the first 
million years or so of disk evolution. This second con- 
straint may be in agreement with Stardust results. Early 
studies found no aqueous minerals in the Stardust mate- 
rials, suggesting that the outward transport and incor- 
poration of these materials into the 81P/Wild 2 come- 
tesimals occurred very early, before the primary accre- 
tion an d fragmentation of the chondritic-asteroid parent 
bodies (|Woodenll2008D . However, the possible identifi- 
cation of igneous materials in the Stardust samples may 
contradict this scenario or else put interesting time-line 
constraints on the formation time-scales of planetesimals 
and cometesimals across the rang e of Solar Syste in radii 
(jStodolna. Jacob fc Lero^ l2010t iJoswiak et all [SoTo). 
Another constraint linked closely with time is the large- 
end grain size for which our model predicts measurable 
outward mixing. Not only do the peak (and midplane- 
flow-upper-limit peak) values of Cn beyond 25 AU drop 
off markedly for grains a few millimeters in size, but the 
time-frame for that outward transport is also restricted 
(to barely more than 10^ years for 2 mm-sized grains), 
so that comet-formation would have to be quite rapid to 
capture such large inner-disk grains if they did appear in 
the comet-forming regions. 

The transport of different-sized particles is also of 
interest for models that include the effects of grain 
growth. Some observational studies of disks suggest 
a link between disk crystallinity and gr ain growth 
(|van Boekel et al.ll2005HOlofsson et al.ll2009[ ). with char- 
acteristic crystalline grain sizes of a few microns (dis- 
tinctly larger than typical ISM grains). At the very least, 
grain growth and crystallization likely occur on similar 
time scales. The results presented here suggest that as 
dust grains grow, they should become less likely to enter 
and more likely to leave the outer regions of the disk. 
Large grains become mostly conflned to small disk radii 
where the disk gas is still dense enough to support them. 
We know that growth at centimeter and meter scales 
poses a constraint to disk models and planetesimal for- 
mation, because particles of this size are expected to fall 
inward toward the parent star very rapidly, depriving the 
disk of these solids. However, growth to millimeter scales 
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Fig. 13. — MaxiinuirL values of normalized concentration beyond 
25 AU as a function of Schmidt number for a range of grain sizes 
and two source regions. Low diffusivity (large Sc) tends to re- 
strict particles from reaching the outer disk. Dashed lines with 
black markers denote the midplane-fiow simulations. Nominal disk 
model. 

may also place constraints on cometesimal formation if 
the outer regions of the disk cannot support particles of 
even these sizes. At the very least, it could restrict the 
composition of bodies formed in the outer disk, if most of 
the particles available are those falling inward from large 
AU. 

5.3. Transport Varying the Diffusivity 

Next, we explore how the outward mixing of inner- 
disk grains depends on the Schmidt number, which is 
the ratio of the disk viscosity to the gaseous diffusivity 
(so that a low er Schmidt number means more rel ative 
diffusivity.) In iPavlvuchenkov fc Dullemondl (|2007[ ) the 
authors argue for an Sc lower-limit of 1/3. Therefore, 
we have run simulations of particle transport varying the 
Schmidt number by a factor of three, and in Figure [T3] 
we plot the values of maximum normalized concentra- 
tion beyond 25 AU for four particle sizes and two source 
regions as a function of Sc. The primary effect shown 
in Figure [13] is that a higher diffusivity can lead to a 
substantially higher degree of outward transport in the 
general accretion-flow case of our particle-transport sim- 
ulations. Here we see variations of an order-of-magnitude 
or more in the peak outer-disk normalized concentrations 
between the highest {Sc = 1/3) and lowest [Sc = 3) dif- 
fusivity simulations across the range of particle sizes con- 
sidered. This is les s than the variation predicted in the 
model scenario of IPavlvuchenkov fc Dullemondl (|2007[ ) 

(crsc=i/3/o'Sc=3 « (R/Rsource) ), but is stih quite a sub- 
stantial effect considering that our normalized concen- 
trations are reported for the whole of the outer disk and 
that our pool of contaminants is limited to those present 
in the source regions at i = only. 

More generally though, the Schmidt number controls 
the relative importance of advection versus diffusion of 
particles. When the diffusivity is low, advection dom- 
inates, so that more grains reach the outer disk if the 
bulk flow is outward, and fewer do so if their bulk flow 
is inward (including the case of 2 mm particles in the 
outward-flowing-gas midplane-flow case whose radial ad- 
vection is dominated by headwind drag). Conversely, 
high diffusivity pushes the system in the direction of a 
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concentration of 20 fim particles beyond 25 AU as a function of 
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est pre-peak time of half-maximum concentration in the accretion- 
flow simulations. Upper curves mark the earliest post-peak time 
of half-maximum concentration for the accretion-flow (solid lines) 
and midplane-fiow (dashed lines) simulations. Nominal disk model. 



flow-independent state where diffusion, both within the 
disk gas distribution and inward onto the parent star, 
is the dominant term. The diffusivity can also have 
an important effect on outer-disk contamination from 
the inner-most source regions. For 20 /Ltm-sized par- 
ticles from the inner-quarter source region, Sc = 1/3 
simulations produce a peak normalized concentration 
in the comet-forming region of 1.7%, up from 0.2% in 
the baseline Sc — 1 simulations. In the low-diffusivity 
simulations inner-quarter source grains of all sizes are 
completely absent from the comet-forming region unless 
directly advected there by outward-flowing gas in the 
midplane-flow simulations. 

Finally, to some degree, the diffusivity also affects the 
timing of outward mixing. This is plotted for 20 ^m-sized 
particles in Figure [T31 The higher diffusivity simulations 
allow higher normalized concentrations in the outer disk 
to occur sooner and to last longer (~ 4 x lO'^-l.l x 10^ 
years half-max to half-max) , whereas for lower diffusivity 
the outer-disk normalized concentration peaks later and 
more briefly (^ 9 x 10^-5.5 x 10^ years). 

Our results for the variation in transport efficiency 
with Schmidt number suggest that two distinctly differ- 
ent disk models could be consistent with observations. 
One possi bility is that the diffusivity i s high , as sug- 
gested bv IPavlvuchenkov fc Dullemondl ()2007D . In this 
case, it seems plausible that a significant fraction of 
grains, formed in high temperature regions close to the 
star, can end up in the comet-forming region, even if 
the disk is highly turbulent and no significant settling 
occurs. The recovery of CAI-type grains by Stardust 
~ whose condensation requires especially high temper- 
atures - poses the strongest constraint on this scenario, 
and we have not demonstrated explicitly that it is pos- 
sible. However, given that the early disk was certainly 
hotter than the disk which we have modeled, it is plausi- 
ble that enough material from the innermost disk could 
reach large e nough radii. Hot outflows near the disk mid- 
plane, which JTscharnuter et al.l ()2009|) find very early in 
the star/disk formation process, would assist outward 
transport from the innermost disk. 

A qualitatively distinct disk model with low diffusivity 
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Fig. 15. — For the static disk models, the normalized concen- 
tration of 20 /im particles beyond 25 AU for two source regions. 
Pale-curves denote the staticO disk model. Dark, heavy curves de- 
note the steady disk model assuming a gas mass truncated beyond 
100 AU. The curves for the nominal, evolving disk model are in- 
cluded for comparison. Simulations run in the accretion-flow case. 

is also conceivable. Our results suggest that if Sc is signif- 
icantly larger than unity, there is negligible transport of 
crystalline silicates subject to the accretion flow to the 
region where Jupiter family comets form. This is true 
even if silicates can form at radii as large as 10 AU. A 
low diffusivity disk, however, could potentially be favor- 
able to the establishment of an outward midplane flow, 
which i s able to move settled particles outward efficiently. 
iCieslal ()2009() showed that a transport scenario domi- 
nated by high-altitude-inward and low- altitude-outward 
advection can lead to segregated grain populations with 
different processing histo ries. This is possibly c ompatible 
with disk observations of lOlofsson et al.l ()2009D that sug- 
gest some disks are more crystalline at larger distances. 

5.4. Transport in More/Less Compact Disks 

Finally, we assess the role of the disk mass distribu- 
tion and its evolution on outward transport. First, we 
contrast the results of runs assuming a static disk with 
those of the evolving disk, and then move on to exploring 
the effects of the initial compactness of the disk in the 
evolving-disk model. Of our two static-disk models, the 
first uses the t — profile of the nominal-disk held static 
(the staticO-disk model). The second is the steady-disk 
limit of the t — nominal-disk profile, where the disk 
surface density follows Equation ([3]), with Rd = oo (the 
steady-disk model). 

Results from the two static-disk models are shown in 
Figure I15[ which plots the normalized concentration of 
20 /im particles beyond 25 AU with time for two differ- 
ent source regions. The mass per AU of the gas disk does 
not drop off in the outer disk for the steady-disk model 
but is instead continuous out to the edge of our simu- 
lation space. Therefore, the normalized concentrations 
beyond 25 AU for this disk model are calculated using 
the disk gas mass only out to 100 AU (for a total 0-100 
AU disk mass of 0.162 Mq — 0.03 Mq of gas is contained 
within the first 21 AU in the steady-disk model). Fig- 
ure [TS] includes the results for the evolving nominal-disk 
model for comparison. The most obvious result from 
Figure [15] is that the two static-disk models, together 
with the nominal evolving-disk model, form a hierarchy 



in outward mixing that is based on disk structure. In 
the staticO-disk model the region of outward-flowing ac- 
cretion flow remains fixed relatively close to the central 
star, and those outward velocities are fairly rapid. There- 
fore, more than twice the peak fraction of nominal-disk 
particles move to the outer disk (and stay there, as per 
the static-disk-model tendencies seen in Section 2]) . The 
surface-density structure in the steady-disk model means 
that in the accretion-flow case the gas radial- velocity is 
everywhere inward, and hence the number of 20 /im par- 
ticles beyond 25 AU peaks at less than 40% of the value 
for the corresponding nominal evolving-disk simulations, 
and those particles eventually all fall back inward onto 
the parent star, — and more rapidly than in the evolv- 
ing disk simulations. Roughly 23% remain in the disk 
after 10^ years (when 30% remain in the nominal-disk 
accretion- flow simulations) . 

However, in terms of the normalized concentration of 
inner-disk particles in the outer disk, the disk-gas mass 
distribution also really sets the steady-disk simulations 
apart from the rest. With so much more mass in the 
outer disk than the other two models, the steady-disk 
model is that much less intrinsically capable of contami- 
nating the outer disk with inner-disk grains than are the 
disks with a more compact disk mass distribution. Fi- 
nally, for the upper limits placed by the midplane-flow 
simulations both static-disk models send more than 90% 
of the simulation particles from all source regions to be- 
yond 25 AU; the outward advection is only more rapid in 
the staticO-disk model. Therefore, the upper bounds on 
the mixing in these models are set only by their relative 
source-region and comet-forming region disk-gas mass ra- 
tios. 

Although some trends are the same for static as com- 
pared to evolving disk models, the results imply that 
the evolution of the disk - particularly at early times - 
has a substantial impact on particle transport. Steady- 
disk models tend to underestimate the degree of out- 
ward particle transport possible in protoplanetary disks 
in the early stages of disk evolution, when the most 
processing of high-temperature materials occurs. The 
outward transport of particles is substantially more ef- 
ficient at early times, consi stent with t h e res ults of 
iBockelee-Morvan et all (|2002D and ICieslal (|2009[) . even 
here assuming a fixed diffusivity. How strong the effects 
of evolution are will depend, of course, on the details of 
the initial conditions for the disk. To examine this depen- 
dence, we vary the initial compactness of our evolving- 
disk model (parameterized by Rd in Equation [3]) . We 
consider _Rd = 5, 10, and 40 AU, in contrast to Rd = 
20 AU in the nominal-disk model. We retain 0.03 Mq 
as the starting mass of the disk, so that the lifetimes of 
the model disks vary between 3.6-6.9 Myr, and the t = 
mass accretion rates between 8.3 x 10~^-8.6 x 10~^Mq 
yr~^. We do not consider possible variations in heating 
for disks forming in more or less compact configurations, 
reta ining the static temperature profile described in Sec- 
tion [SUEl. 

^ We also continue to employ the same uniform radial distribu- 
tion to initiate our simulation particles. One should note, however, 
that the gas-mass-per-AU is no longer roughly uniform across the 
source regions in the ij^j = 10 and 5 AU disks. For Rj = 10 AU, the 
ratio of the average-gas-mass-per-AU for the inner-quarter source 
region to the outer-half source region is ~ 1.5, and for _R^ = 5 AU 
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Fig. 16. — Maximum values for the normalized concentration of 
20 iim particles beyond 25 AU as a function of initial disk compact- 
ness. Dashed lines with black markers denote the midplanc-flow 
simulations. Grey line marks the baseline case of Rj = 20 AU. 
Note that the most compact disks both send a larger fraction of 
simulation particles beyond 25 AU and have a larger fraction of 
disk material in the source regions at i = 0. 
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Fig. 17. — Maximum values for the normalized concentration of 
2mm particles beyond 25 AU as a function of initial disk compact- 
ness. Dashed lines with black markers denote the midplane-flow 
simulations. In the accretion-flow case, almost no 2mm particles 
from the inner-quarter source region reach beyond 25 AU. 



Reducing the compactness of the model disk has two 
consequences: the disk initially expands more rapidly, 
and it does so from an initial state where a greater frac- 
tion of the mass is at small radii (and thus hot and poten- 
tially able to form crystalline material). The combined 
magnitude of these effects is shown in Figure [TBI where 
we plot the maximum normalized concentration of 20 fim 
particles beyond 25 AU as a function of Rd- Variations 
of approximately one order of magnitude in R^ result 
in nearly two orders of magnitude difference in the peak 
concentrations for each source region. For the smallest 
disk - with i?d = 5 AU - the normalized concentration 
of all simulation particles peaks at ~68%, more than an 
order of magnitude in excess of the value for the nom- 
inal simulations. A substantial part of this increase is 
due to the larger reservoir of mass that a compact disk 
has at small radii, which can potentially contaminate 
the (proportionally less massive) outer disk. That this 
mass-distribution effect is important is made clear by 
the fact that even the upper-limit concentrations from 
the midplane-flow simulations show a strong (order-of- 
magnitude) trend with Rd- 

The fact that the outer disks of the initially more- 
compact-disk models are less massive (and more tenu- 
ous) does place some constraints on the outward mixing 
of large grains in our particle-transport simulations. In 
Figure [17] we plot the peak normalized concentration of 
2 mm particles beyond 25 AU as a function of Rd- All 
of the disks, of course, have less outward transport of 
these larger grains, but in our simulations, the peak in 
the fraction of particles sent to beyond 25 AU occurs 
for Rd = 10 AU, which sends up to 0.25% of all 2 mm- 
simulation particles beyond 25 AU (4% in the midplane- 
flow case) compared to less than 0.1% for Rd — 5 AU 
(less than 2% in the midplane-flow case). Because of 
the larger source region mass for Rd = 5 AU, this drop 
in fractional outward transport constitutes a leveling off 

it is ^ 2.6. This afli'ects the accuracy of the normalized concentra- 
tion values that we report for these disks (especially for the widest 
source regions), but it does not alter the trends with ij^ repeated 
across source-region populations. 
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Fig. 18. — Mass accretion rates for disk models initiated with the 
same total disk mass, but with different initial exponential cutoff 
radii (Rd)- Note the reversal in ordering from early to late times. 
All disks were initiated at t = with a total mass of 0.03 Mq. 

in the maximum-achievable normalized concentrations of 
inner-disk particles in the outer disk. 

Comparing these results to observations of crystalline 
silicates in other disks allows us to construct a possi- 
ble scenario for the formation and transpo rt of high- 
tempe rature minerals in those disks. In the W atson et al.l 
(|2Q09D survey, the disks observed were around T Tauri 
stars ~l-2 Myr old and showed a slight correlation be- 
tween the crystalline-silicate mass fraction and the mea- 
sured accretion rate onto the star. As shown in Figure 
[THJfor our model disks, t — 1 Myr corresponds to a time 
when those disks that were initially most compact (with 
high t — accretion rates) now have lower accretion rates 
than the initially more-extended disks. Therefore, small 
grains that formed at small AU near t = should be 
more broadly distributed in disks that now have lower ac- 
cretion rates. However, lower accretion rates also mean 
disks that are more tenuous. The initially most com- 
pact disks lose mass the fastest, causing their outward 
transport efficiency to drop rapidly with time. There- 
fore, grains that formed at small AU after t ^ 2x 10^ yr 
may be better distributed in the initially more-extended 
disks, which now have higher accretion rates. 
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These properties could potentially explain the oppos- 
ing trends betwe en crystallinity and o bserved accretion 
rate reported by iWatson et alj (|2009[ ). Pyroxene is fa- 
vored over olivine when the minerals are formed by 
condensation (IGaill 120041 : IWooden. Harker. fc Brearlevl 
[2005.) . Therefore, a possible scenario is that pyroxene 
formed primarily at very early times when the disks were 
hottest and became most broadly distributed throughout 
the initially most-compact disks. Most of the olivine, on 
the other hand, may have formed over a longer period 
of time. It might, therefore, have become more broadly 
distributed within disks that retained significant surface 
densities for longer times (the initially more-extended 
dis ks). This hypothesis does not, however, explain why 
the IWatson et al.l ()2009[ ) survey found no correlation be- 
tween crystalline mass fraction and disk mass. 

6. CONCLUSIONS 



Observations of comets (iHanner et al.l 119941 : 
IWooden. Woodward, fc Harked I2004D and sample 



return from the Stardust mission (jBrownlee et al.ll2006D 
suggest that outward transport of high-temperature 
solids to the comet-forming regions is a necessary 
component of any successful model of the early Solar 
nebula. We have constructed a ID model of viscous 
disk evolution and particle transport that we use 
to study the relationships between the evolving disk 
mass distribution and the local conditions that are 
supportive of the outward mixing of grains, specifically 
including self-consistent treatment of the dust grains' 
surface-area-to-mass ratios in our transport terms. We 
have examined patterns of outward mixing considering 
variations in the radial gas-flow structure, the sizes of 
dust grains, the model's Schmidt number, and the initial 
compactness of the model disk's gas mass distribution. 

We find that a range of disk models are able to 
account for the presence of high-temperature-inner- 
disk particles in the com e t-form ing region. As in 
iDullemond. Apai. fc WalchI (|2006D . the most favorable 
models involve rapid early expansion of initially quite 
compact disks. These conditions can result in inner- 
disk material flowing to th e outer disk o n a short time 
scale. In agreement with iCieslal ()2009() . we find that 
outward transport is more efficient at early times due 
largely to stronger gas outflows. Even more favorable 
are cases involving outward-flowing gas at the disk mid- 
plane (Keller & Gail 2004; Cicsla 2007, 2009), which are 
capable of delivering substantial fractions of grains in- 
ward of 1 AU to the outer regions of the disk. It is not 
known whether such outward flows exist in real disks (or, 
whether the particles are sufficiently settled to experience 
them consistently) . Whether such midplane fiows are re- 
quired to explain the Stardust results depends upon the 
Schmidt number, which at its lower limit (highest rela- 
tive diffusivity) could allow for a relatively high degree 
of general outward mixing, even from relatively near the 
parent star. While the mean accretion flow is less effi- 
cient at transporting particles outward, significant out- 
ward flows of crystalline material would be possible if 
the disk were initially both compact and able to form 
silicates out to a distance of several AU. 

Our results suggest that it is the largest particles re- 
covered by Stardust that place the tightest constraints on 
the evolution of the early Solar nebula. Particles as large 



2 /zm are well coupled to the disk gas, and plausible lev- 
els of turbulent diffusivity allow them to reach the outer 
disk from small radii in a wide variety of models. Our 
simulations show that 20 //m particles can also reach the 
outer disk — consistent with the Stardust results — as 
long as they were incorporated into larger bodies within 
the flrst 1-2 Myr of the disk evolution. For even larger 
particles, those of a few mm or larger, none of our disk 
models admit significant outward transport to the comet- 
forming region, regardless of the proposed source region. 
Even grains of a few hundreds of microns in size should 
be rare, since their residence time in the outer disk is 
very limited. Discovery of such grains would, within the 
context of the model developed here, require a more mas- 
sive disk during the time of their formation, transport, 
and incorporation into cometesimals. 

The existence of substantial uncertainties in the disk 
physics currently precludes an observational determina- 
tion of the characteristic radii at which crystalline silicate 
formation occurs. The strongest constraints come from 
observations of CAI-like grains, which probably formed 
at fractions of an AU from the Sun. Only a subset of 
our accretion-flow simulations arc likely compatible with 
this finding. A broader range of models in which the mid- 
plane gas flow is outward are viable, since these models 
retain particles in the disk for a longer period of time 
and facilitate the co-existence of grain populations of 
markedly different processing histories. Generically, we 
expect that provided sufficient outward transport is pos- 
sible, turbulent diffusion will cause particles to become 
radially well- mixed over tens of AU so that incorporation 
of CAI-like grains and Fe-rich crystalline silicates (possi- 
bly formed in water-rich shocks out to 10 AU) into the 
same body is plausible. 

For astronomical observations, our results suggest that 
the mass and size of the disk at an early epoch are crit- 
ical parameters in determining the viability and extent 
of outward particle transport. These disk properties are 
inherited from the mass and angular-momentum struc- 
ture of the cloud that collapses to form the star, and 
would vary from system to system. We expect this to 
result in variations in the crystalline fraction at large 
radii in disks around different stars. We have also noted 
that the time scale during which conditions allow sub- 
stantial outward transport is only a small fraction of the 
disk lifetime. A limited time for outward mixing is both 
interesting and puzzling when paired with the emerg- 
ing evidence for aqu eous and/or igneous grains among 
the Stardust samples (lWoodenll2008l: iJoswiak et al.ll2010l: 
iStodolna. Jacob fc Lerouxll2010D . However, it may allow 
the Stardust results to be reconciled with other obser- 
vations that require compositional gradients within the 
disk to be maintained. A possible scenario is that an 
early stage of rapid expansion drove particles into the 
outer disk where they were rapidly locked up into larger 
icy bodies. Once the disk had expanded, the gas flow 
became more uniform, limiting the outward transport of 
solids. 
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APPENDIX 
VERIFYING THE FIDELITY OF THE PARTICLE TRANSPORT MODEL 

Our particle-transport model consists of two pieces. The first is the simple mean-radial velocity of particles immersed 
in a radial and azimuthal gas fiow. It is represented by a calculation of the mean-radial velocity at each grid point in 
the model disk at each time step. The second piece is the turbulent diffusion of the particle ensemble. It is represented 
as a random walk with respect to the background mean-radial velocity. For the random walk, we calculate a turbulent 
velocity for each particle at each time step and the probability that the particle will randomly walk inward versus 
outward. Using these velocities and random-walk probabilities, we can then do a simple integration of the radial 
motion of each particle using Equation ([T]) . 

In this appendix, we demonstrate that this simplified approach to particle mixing and transport does, in fact, 
accurately represent the physics necessary to our particle-transport model. In Section lA.li we demonstrate that we 
have chosen appropriate local velocities for the dust-mean-radial velocities and that those velocities are sufficient 
to represent radial particle motion. We compare our radial-dust trajectories t o tr ajector ies produced usi ng di rect 
numerical integration of the the dust-particle equations of motion, Equations (|A1[) and (jA2[) . In Section [A. 21 we 
compare test-case simulations from our particle- transport model to analytical solutions of the diffusion of a contaminate 
in an accreting disk given by iClarke fc Pringlj (fTQSS.l . We demonstrate that the collective behavior of our random- 
walking particles produces the equivalent of diffusion in a gas disk. 

Advection Compared to Force-Balance Integration 

Our particle-transport simulations use an approximation technique to calculate particle trajectories. Here we com- 
pare the trajectories produced using this technique to precision numerical integrations of those trajectories. T he 
equations of motions for a particle in a laminar disk-gas flow in the Epstein-drag regime are (|Takeuchi fc Linll2002f) 



"r.d _ "Id - "k 



-Cr Pgl^thcrm iVr,d " ^^r,g) , (Al) 



dt R 3 

dv,p,d «0,d«r,d 1^ , ^ ..„s 

—jp = ^ -^Cr PgWthorm (V4,A " W0,g) , (A2) 

where Ur,d and u^^d are the radial and azimuthal velocities of the dust particle; Vr,g and v^^g are the local radial and 
azimuthal velocities of the disk gas; vk is the local Keplerian velocity; Cr is the particle surface-area-to-mass ratio; 
and pg and Wthcrm are the local-gas volume density and thermal velocity, respectively. 

In the case of laminar-gas flow with no turbulent forcing of the particle motions, we can numeri cally solve for each 
particle trajectory directly. To do this, we use a Runge-Kutta-style integrator, where Equations (|Al|) and (jA2|) are 
coupled with 

—7- = Vr.d (A3) 

at 

to define the particle motion radially within the disk. The integrator uses step-size adjustment to reach a relative 
accuracy of 10~^-10~^ in rj, Wr,d, and w^.d. Note, however, that the integrator becomes numerically expensive when 
the individual terms in Equations (JA1[) and (|A2p are large but sum to values that are small. Therefore, if the initially 
computed acceleration {dvr.d/dt or dv^^^/dt or both) is very small (< 10^^"* cm s~^), the integrator sets the small 
acceleration equal to zero. Also, if w^.d ^ vk (10"^*^ relative difference) then the vi ^ — w^-term in Equation (jAip is set 

equal to zero, but only if -C/jpgUthorm > -pr- If W0,d ~ wk but the second condition is not met, then the wl ^ — w^"term 

3 i? ' 

is still important and the integrator instead calculates the particle trajectory to the lower relative accuracy of 10~^. 

Most of the comparison simulations presented below are run using the midplane-flow-velocity case, as defined in 
Section 13.21 of this paper. We begin by presenting particle trajectories calculated within an initially tenuous disk 
(Rd = 500 AU, Sg (lAU) ~ 1 g cm~^ at t = 0) that is photoevaporatively cleared within <10^ years. Radi al pa rticle 
trajectories in such disks can be complex and so provide good tests of our particle-transport model. Figure lAlJ plots, 
for both numeric integrations and particle-transport simulations, the trajectories and radial velocities of 0.2 mm-sized 
particles initiated at three different radial locations within this tenuous disk. Here, such large particles are only 
loosely coupled to the gas flow, with exponential-stopping times (T^) at i = of a few to greater than ten Kepler 
times. Nevertheless, there is good agreement between our particle-transport simulations and the more precise numeric 
integrations. 

However, our particle-transport simulations do lack certain details; as shown in Figure IA2| they do not account 
for orbital eccentricity. Th us, the particle-transport simulations are limited in their depiction of complex gas-particle 
interactions. In Figure 1X2] these interactions pertain to a particle disengaging from and outward-sweeping photoevap- 
oration front. However, the results of the particle- transport simulations presented in this paper are not sensitive to 
such a high level of detail. Instead, we are most interested in the broad effects of disk mass and particle size. For the 
tenuous disk, we ran simulations for particles of sizes 0.2 /xm, 2 /zm, 20 /xm, and 0.2 mm from 20 starting locations, 
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Fig. A1. — Radial position (a) and radial velocity (b) of simulated particle trajectories. Comparison between precision numerical 
integrations (heavy, dark lines) and the particle-transport-model trajectories (thin, pale lines). In this figure, 2 mm-sized particles were 
initiated at three different radial positions within a tenuous, evolving disk and were subject to the midplane-flow gas-velocity case. 




m 


1500 


U 


1000 


>i 


500 


-P 

-H 

o 
o 





0) 

> 


-500 


rH 


-1000 


-H 





5x10 6x10 7x10" 8x10' 9x10' 
Time (years) 



(^ 5x10 




6x10 



7x10 



8x10" 



9x10 



Time (years) 



Fig. A2. — Radial position (a) and radial velocity (b) of a 2 mm particle initiated at 242 AU within a tenuous, evolving disk in the 
midplane-flow gas-velocity case. Comparison between a precision numerical integration (black) and the particle-transport model trajectory 
(red). 

0.5-500 AU. The exponential-stopping times in the simulations ranged from less than 10^* to greater than 10 Kepler 
times, and in general, the relative radial positions of the simulated trajectories agreed with the numeric integrations 
to within 2% or better. 

Where our particle-transport simulations deviate the most from the numeric integrations is for particles experiencing 
a long period of steady infall. This is true regardless of how well coupled the particles are to the gas motions, and 
is an effect of the particle-transport model repeatedly overshooting the estimate of the infalling particle trajectory. 
Figure IXsT a) plots trajectory comparisons for particles initiated at 500 AU in a massive, but initially compact disk 
(Md = O.O5M0 and Rd — 10 AU). In this scenario, the particles begin very decoupled from the gas (Tg > 10"^ Kepler 
times), and after their initial infall toward the main disk they become moderately well coupled (Tg ~ 10~^ and smaller). 
In Figure lASf b) we plot the relative difference in the particle radial positions between the particle-transport model 
trajectories and the numerical integrations. Here we can see that for each infall, these two trajectory calculations 
diverge as the simulated trajectories outpace those of the numeric integrations. Therefore, exaggerated, rapid infall of 
particles may be a drawback of our particle-transport model. However, inspection of all of our simulated trajectories 
suggests that this exaggerated infall should have a negligible effect on our general understanding of the timing and 
trends of particle transport. As a check we have also run comparison integrations using the accretion- flow case in our 
nominal-disk model (Section l5.1l) . We ran four particle sizes (0.2 /im-0.2 mm) initiated at 15 AU and found positional 
differences of less than 2% for all but the very final stages of infall. 

Therefore, in the absence of turbulence, the simulated trajectories produced by our particle-transport model appear 
to be a sufficiently accurate representation of radial particle transport in an evolving disk under the influence of gas 
drag, both from the pressure-supported-azimuthal and the radial gas flow. 

Turbulent Diffusion Compared to Analytical Solutions 

Here we demonstrate the fidelity of the random- walk method in our particle-transport model in reproducing accurate 
diffusion of the particle ensemble. Using the analytical solutions of .Clarke fc Pringlc, (|1988. 1 we compare the time- 
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(a) Particle trajectories initiated at 500 AU within a massive, 
midplane-flow gas-velocity case. Comparison between precision numerical integrations (heavy, dark lines) and particle-transport model 
trajectories (light lines), (b) Relative difference over time in the radial positions between the numerical integrations and the particle- 
transport simulations, [(simulated - integrated) / integrated] Note that the time-axes are in log-scale. 

dependent particle distributions from our particle-transport simulations to the analytically expected contaminant 
distributions and discuss the effects of our choice of simulation time step on the results. All solutions and simulations 
pres ented below use stetic-st eady-disk-surface-density profiles. 

In iClarke fc Pringlj ()1988l ) , the authors solve for the time-dependent concentration, C, of a contaminant initiated 
in a ring at specified radius: C (i = 0) = CqS {R — Rq), where C is the ratio of the local contaminant mass over the 
local disk mass, R is the distance from the central star, and Rq is the initial position of the contaminant. The disk 
that they consider has a surface-density profile given by E = Egi?"", where Eg and a are constants; it is also a steady 
disk, so that for disk viscosity, i' = voR^ {b and i/q constants), b = a, and the gas-accretion velocity may be written as 
vr = voR''~^ where vq = —iua/2. 

Assuming azimuthal symmetry an d that the gas and contaminant are vertically well mixed, the diffusion equation 
in polar coordinates may be written ()Gail|[200ll ) 



d Id Id 



RDig E 



dC_ 
dR 



P,{R,t) 



(A4) 



where Dig is the coefficient of diffusion of the contaminant within the gas, an d Pi corresponds to the rate of production 
of the contaminant per radial increment within the disk. IClarke fc Pringlg ^1988) (and the setup of our simulations) 
assume zero production of the contaminant, Pi{R,t) = 0, and that the diffusivity scales with the disk viscosity. 



Dig — C^i where C is a constant. For these particle-transport simulations using a static disk profile. Equation (jA4p is 
further constrained by dT,/dt = and so vr refers specifically to the set dust-mean-radial velocity. Using the profiles 
for the disk surface density and viscosity above, the full diffusion equation applicable to these simulations is: 
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which simplifies to 

J_d£ 

vo dt uo 

in the steady-disk (a — b) scenario of IClarke fc: Pringli (|1988D . The scenario of IClarke fc Pringli (|1988[) assumes an 
evolving (accreting) gas disk, but, because it is steady (has an unchanging relative-surf ace-d ensity profile) and the 
diffusion equations do not depend on the actual value of the gas surface density. Equation (|A6|) fits both that scenario 

and our static steady-disk scen ario. 

The general solution given in IClarke fc: Pringlg (|1988l ) to Equation (jA6P using the initial contaminant distribution 
given above is 
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where Ip {x) is the modified Bessel function of the first kind of order 
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However, Equation (|A7I) is not valid for the case of a == 6 = 2, for which [Clarke fc Pringld (|1988f) provide the solution 
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(A8) 



Note, however, that Equation (jASp does not match Equation (3.1.3) in IClarke fc Pringli (J1988I) because of a typo 
in that paper, which accidentally presents the solution for twice the background velocity, v^ — > "^vr- Also, when 
a = b = 2, the natural space to solve the equation is logarithmic space. This happens to be the space used for the 
width of our model-disk grid cells, and may be related to the only odd behavior that we find in our particle transport 
simulations (see below). 

From our particle-transport simulations, we output the number of particles per radial grid space. Therefore, to com- 
pare the solutions provided by Equations (jA7P and (jASP to our simulations, we must first convert the analytic C {R, t) 
into an expected fractional-mass distribution. The concentration can be expressed as the ratio of the contaminant-to- 
gas surface densities, C — crfS, where a is the surface density of the contaminant. Because these test cases assume a 
static-disk-surface-density profile, it is then relatively simple to solve for the fractional mass of the contaminant rrii, in 
a given grid space, where i denotes the grid space, which we will normalize by mtot, the total mass of the contaminant 
in the disk at t = 0, where 
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= [RoT. {Ro) Co] ^ I i?S (i?) C {R, t) dR , 

ymtotj Jr=r.-_-^/^ 

where Ri-1/2 and Ri+1/2 are the inner and outer boundaries of the grid space, respectively. Combining Equations (|A7p 

and (|A9[) . and using a variable substitution, s — i?^'', we can write the general expression for the expected mass 
distribution as 
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We solve the integral numerically using a Simpson integrator and the Numerical Recipes functions for the modified 

Bessel functions (Press ct al. 1992). 

For the case of a = 6 = 2 given in Equation (jASp the expected mass distribution in a static disk is given by 
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where Y — ln(i?/i?o), and we again, we use a Simpson integrator to solve for the expected {mj/ mtot) for t hese c ases. 

To compare our particle-transport model to the analytic test cases described in Equations (|A10p and (|A11I) . we 
ran simulations of 10,000 particles initiated at 40 AU; we include some simulations with zero background velocity, for 
which we substitute vq — into the solution equations above (instead of the accretion-case vq = — 3i'o/2). We vary 
the value of a, but all simulations use v (lAU) = 4.941 x 10^^ cm^ s~^ and C = 1- In these test-case simulations, the 
values of Eg and the par ticle size are arbitrary, since the background velocity of the particles, vr, is prescribed. 

Panel (a) of Figure lA4l plots a comparison between our numerical simulations of particle transport and the analytically 
expected contaminant-mass distribution for the case of a = 1.5 and zero background velocity at two different times. 
It shows that the overall, evolving distribution of the particle ensemble is in good agreement with the analytical 
solution. However, the scatter produced by the discreet nature of the particle simulations makes it difficult to judge 
the accuracy with which the simulations reproduce this analytical test case. Therefore, Figure IA4I panel (b) plots 
an alternate comparison of the simulated and analytical mass distributions. Here the values plotted correspond to 
a summing of the fractional-mass per grid space either from R = outward, or from R = Rq = 40 AU inward (for 
R < Rq) and outward (for R > Rq). Again, the simulations and the analytical solutions are in good agreement. 

Simulating particle diffusion with zero background velocity in a steady disk for a range of a-values, we find, in 
general, that our simulations agree well with the analytical solutions. The exceptions occur for very steep (a = 3), or 
very shallow (a = 1) disk profiles and are caused by the finite nature of our simulation space. A viscosity (diffusivity) 
profile that goes as R^ leads to rapid transport of particles over very large distances at large R. Therefore, the outer 
edge of our simulation space acts as a sink for particles. At early times, before a significant fraction of particles have 
been lost, the simulations and analytical solutions are in good agreement, as shown in Figure Esl they deviate from 
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Fig. A4. — (a) Comparison of the mass distributions of the particles (contaminant) between the particle-transport simulations (pale lines) 
and the analytical solution (dark lines) at two times. The vertical axis is the fraction of total particles (total contaminant mass) per radial 
grid space, (b) Comparison between simulations (pale lines) and analytical solutions (dark lines) of the summed-mass distributions. Solid 
lines represent the distributions summed from /? = 0, and dashed lines represent the distributions summed from K = Ro = 40 AU. For this 
figure a = b = 1.5 and no = 0. For the simulations, 10,000 particles were initiated at 40 AU and diffused using a global time-step size of 
At = 2000 years. 
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Fig. A5. — Mass-distribution profiles (a) and summed profiles (b) at three times for the case of a = 6 = 3 with vq = 0. Comparison 
between simulations (pale lines) and analytical solutions (dark lines). In (b) solid lines represent summing from R = 0, dashed lines from 
R = Rq. In this simulation at i = 100 years, 100% of the particles remain on the grid, but at t = 500 and 2000 years only 93% and 87% 
remain, respectively. The simulation was run using a global time-step size of At = 100 years. 

each other at later times. In the a = 1 case, it is loss past the inner boundary of the simulation space that causes the 
simulations to deviate from the analytical solutions. 

Loss of simulation particles past the simulation-space boundaries is not an issue when we include the background- 
accretion velocity. In Figure EBl we plot the mass-distribution profiles at f = 10^ years for a range of a- values with the 
accretion velocity included. Because vr is inward, it keeps particles away from the outer boundary, and loss past the 
inner boundary is a part of the analytical solution as well as the simulations. All of these cases show good agreement 
between the simulations and the analytical solutions. 

Finally, we would like to be sure that our particle-transport model correctly simulates diffusion of the particle 
ensemble for the global time steps used in our static-disk models (At = 2000 years), and for the time steps required 
by the disk-evolution model (At ~ 0.05 years). As was described in Section [5751 maximum limits are placed on the 
particle time stepping to ensure that disk conditions local to each particle are obeyed. However, smaller time steps 
mean more steps to simulate a given amount of time in the disk; therefore, the probabilities of a particle stepping 
inward versus outward are closer and closer to equal. Cumulatively, however, all time-step sizes used must evolve 
the distribution of the particle ensemble equally. Therefore, we have run a few test-case simulations using time-step 
sizes ranging from 10^"^ years to 100 yea rs. For the case of a = 1 and vq — 0, the resultant distributions versus the 
analytical solutions are plotted in Figure IA7I All of these simulations agree well with the analytical solution, and the 
scatter produced for simulations with different time-step sizes are all comparable. 

While the agreement between simulations using different time-step sizes is good for almost all of the simulations 
we tested, one set of simulations shows a marked change in the scatter depending on the size of the time step used. 
Figure ESJ plots simulations run for a range of time-step sizes for the case of a = 6 = 2 and vq = — 3z/o/2 compared to the 
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Fig. A7. — Mass-distribution profiles (a) and summed-profiles (b) at t = 2 X 10^ years for the steady-disk case where a = b = I with 
VQ = 0. Simulations were run for a range of global time-step sizes (colored lines) and are compared to the analytical solution (black/grey 
lines). In (b) the dark lines represent summing from R = 0, and the pale lines summing from R = Rq. The colored lines of the simulations 
are largely obscured by the analytical solutions (black/grey) plotted over top. 

analytical solutions at two times. Looking at the sumnied-mass distributions, it is clear that all of these simulations do 
distribute the particle ensemble appropriately for this diffusion test case. However, the mass-fraction-per-bin profiles 
show significant, resonant-like scatter in the At = 10 years and At = 1 year simulations; the scatter also appears a 
bit higher than usual in the At — 10^^ year simulation. While it is not clear what causes this large scatter for these 
time steps, it is suggestive that the natural space for diffusion in this regime (a = 6 = 2 — > In (i?)-space) is the same 
as the space wherein our disk-model grid cells are equal width. Furthermore, the resonant-like-scatter region in the 
At = 10 years simulation is restricted to outward of about 15 AU, and this region moves outward with time. Inward of 
15 AU in these simulations, the Kepler time (l/i7K) is less than 10 years and decreases toward smaller R; the time-step 
sizes of individual particles at small AU are restricted to smaller than the global 10-year time step. Therefore, the 
resonant-like-scatter behavior is only observable when the same local time-step size is used for the entire region and 
when the source of all the particles is a single point in R. Nonresonant-like scatter at small AU allows the particle 
distribution to spread out randomly, following the local diffusivity and removing the single-source-point condition 
starting at small AU and moving outward. 

While this resonant-like-scatter behavior is curious and somewhat vexing, it is not a problem for the simulations 
presented in this paper. Aside from the fact that our model disks (Section 13. ip use i' = vqR {h = 1), the time- 
step sizes for our steady-disk and staticO-disk models are large enough that the local time-step restrictions of the 
particle-transport model apply almost everywhere. The evolving-disk models use a time-step size that is small enough 
(At ~ 0.05 years) that resonant-like scatter is not expected to be an issue. Note, though, that resonant-like-scatter 
behavior may be something to watch out for when using this method of particle-transport modeling. 
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